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Abstract 

The aim of the paper is to discuss the main characteristics of a complete theoretical and numerical 
model for turbulent polydispersed two-phase flows, pointing out some specific issues. The theoretical 
details of the model have already been presented [Minier and Peirano, Physics Reports, Vol. 352/1-3, 
2001 ]. Consequently, the present work is mainly focused on complementary aspects, that are often 
overlooked and that require particular attention. In particular, the following points are analysed : the 
necessity to add an extra term in the equation for the velocity of the fluid seen in the case of two- 
way coupling, the theoretical and numerical evaluations of particle averages and the fulfilment of the 
particle mass-continuity constraint. The theoretical model is developed within the PDF formalism. The 
important-physical choice of the state vector variables is first discussed and the model is then expressed as 
a stochastic differential equation (SDE) written in continuous time (Langevin equations) for the velocity 
of the fluid seen. The interests and limitations of Langevin equations, compared to the single-phase case, 
are reviewed. From the numerical point of view, the model corresponds to an hybrid Eulerian/Lagrangian 
approach where the fluid and particle phases are simulated by different methods. Important aspects of the 
Monte Carlo particle/mesh numerical method are emphasised. Finally, the complete model is validated 
and its performance is assessed by simulating a bluff-body case with an important recirculation zone and 
in which two-way coupling is noticeable. 



1 Introduction 



Dispersed two-phase flows, where a continuous phase (a gas or a liquid) carries discrete particles (solid 
particles, droplets, bubbles, . . . ), are of great interest in environmental studies and engineering applications, 
such as dispersion of small particles in the atmosphere or combustion of fuel droplets in a car engine. 

To simulate these flows, the basic equations must be written: the Navier-Stokes equations for the fluid 
phase and the momentum equation for a single particle embedded in a turbulent flow, the latter issue still 
being a subject of current research. For small particle-based Reynolds numbers Rep (whose definition is 
specified below) and particle diameters that are of the same order of magnitude as the Kolmogorov length 
scale, a general form of the particle momentum equation has been proposed HUH. 

In the present work, only heavy particles {pp ^ pf) are under consideration and the equations of 
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motion for a particle can be written: 



^ = "- (la) 

where Us = U(xp(i),t) is the fluid velocity seen, i.e. the fluid velocity sampled along the particle 

trajectory Xp(t), not to be confused with the fluid velocity U/ = U(x/(i), t) denoted with the subscript 
/. The particle relaxation time is defined as 



Pp 'idp 

Pf 3CD|Ur 



(2) 



where the local instantaneous relative velocity is Ur = Up — Us and the drag coefficient Cd is a non-linear 
function of the particle-based Reynolds number, Rcp — dpjUrl/i'/, which means that Cd is a complicated 
function of dp, the particle diameter ||3]- A very often retained empirical form for the drag coefficient is 

{24 
[1 + 0.15i?e°-^s^] ifi?e„< 1000, 
Rep^ p J P _ , 
0.44 if Rep > 1000. 

In the present work, attention is focused on some aspects of the problem. In particular, only dilute 
incompressible gas-particle flows are considered, so that particle-particle interactions are neglected but 
two-way coupling is retained, which means that particle dispersion and modulation of turbulence by the 
particles are accounted for 

The complete problem is formed by the discrete particle equations given above, Eqs ([T]) to (|3]l and the 
field equations of the fluid phase, the continuity and the Navier-Stokes equations, supplemented with a 
source term S that represents the force exerted by the particles on the fluid 



dU 



/J 



= 0, (4a) 



^ + U = + + (4b) 

dt dxj Pf dxi dx"^- 

An "exact approach" (in the spirit of DNS) is possible |4 1, but in practice, the exact equations of motion are 
not of great help. Indeed, in the case of a large number of particles and of turbulent flows at high Reynolds 
numbers, the number of degrees of freedom is huge and one has to resort to a contracted probabilistic 
description. 

Following the classical approach used in single-phase turbulence, one can think of writing directly 
mean-field equations for a limited number of particle statistics (mean velocity, kinetic energy, ...) as in 
fc — e or Rij — e modelling. This is the basis of the Eulerian approach |5 6|. However, due to the complex 
dependence of Tp on particle diameters and on fluid and particle instantaneous velocities, the drag term 
represents a non-linear but local (in term of particle variables) source term. The resulting closure problem 
that appears in the Eulerian approach is therefore difficult. Actually, this issue is very similar to the one 
appearing in the modelling of single-phase turbulent reactive flows |7| and, in this case, PDF models that 
can treat the reactive source terms without approximation have shown their great potential. For the same 
reason, a PDF approach to polydispersed turbulent two-phase flows is interesting. In practice, mean-field 
equations {Rij — e) are used for the fluid whereas a particle pdf equation is solved by a Monte Carlo method 
using a trajectory point of view (Eulerian/Lagrangian models). The PDF model is therefore formulated as 
a particle stochastic Lagrangian model (a set of SDEs). 

Numerous Eulerian/Lagrangian two-phase flow models have been proposed (most of the time with 
interesting and clear ideas), but often with a discrete formulation (in time) and without making the connec- 
tion with a PDF model. When Eulerian/Eulerian {i.e. both phases are described with mean-field equations) 
and Eulerian/Lagrangian models are compared directly, the PDF framework is helpful to reveal that these 
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methods do not contain the same level of information: Lagrangian models are PDF models from which Eu- 
lerian models can be extracted in a consistent way ||6l[8J. The specificity of the present work is to present 
a Lagrangian model, based on a Langevin equation for the velocity of the fluid seen, in the PDF context. 
The theoretical formulation of the Langevin PDF model has already been developed |l8l|9l and the purpose 
of the present paper is to give an overview of the complete theoretical and numerical issues insisting on 
complementary points. In the same PDF framework, an alternative formulation consists in writing a PDF 
equation in analogy with Boltzmann kinetic equation ITOl . that is only for particle location and velocity, 
without considering directly the velocity of the fluid seen. For a comprehensive review of general results 
and methods in particle dispersion, we refer also to Stock's paper 1, 1 1 J . 
More specifically, the aims of the present paper are: 

(i) to outline the main aspects of a PDF model, the interests and limitations of current state-of-the-art 
Langevin models and the keys points of numerical algorithms, 

(ii) to point-out and to address new specific issues for Lagrangian models, such as the addition of extra 
terms for two-way coupling, numerical averages and the mean-continuity constraint, 

(iii) to validate the complete model and to show how it performs by comparing the numerical results with 
experimental ones in a practical case. 

The paper is organised as follows: in Section |2] several mathematical notions related to stochastic 
modelling are clarified, that is the equivalence between the trajectory and pdf points of view, and the 
modelling strategy which is adopted in the present work (the particle-tracking approach). The dimension 
of the system, that is the dimension of the state vector, is also given based on physical principals. In 
Section|3] closure proposals are put forward for the fluid velocity seen, in the form of Langevin equations. 
Emphasis is put on the terms to be added in order to model cases with two-way coupling. In Section |4l 
the numerical approach is presented. The main steps of the particle-mesh algorithm are explained, while 
particular attention is devoted to the problems of defining averages in two-phase flows and of verifying 
the particle mean-continuity constraint. These models are validated in the simulation of a practical case of 
gas-solid flows, Section|5] 



2 Stochastic modelling 
2.1 Mathematical background 

In this section, basic results, concerning the mathematical background of the approach and the correspon- 
dence between SDEs and Fokker-Planck equations, are recalled ll7l [T2l[T3l . If one considers a system of N 
particles interacting through forces that can be expressed as functions of variables attached to each particle 
(for example, position, velocity, ...), then all available information is contained in the state vector, Z, of the 
complete system 

Z_ / /7I /7I rrl , ry2 ryl ryl . . rr N rrN r^W N ( C\ 

— 1^1 : ^2 ' • • ■ ' I ^1 ' ^2 ' • ■ • I ' ■ ■ • I '^l ' ^2 ! • • ■ I J' 

where represents the j-variable attached to the particle labelled i. The dimension of the state vector is 
then d ~ N X p where N is the total number of particles and p the number of variables attached to each 
particle. The complete system, that is the N-particles, is closed. In classical mechanics, the time evolution 
of such systems is often described by a set of ordinary differential equations 

^ = A(i,Z), (6) 
which corresponds, in sample space, to the Liouville equation lfT2]| 

M^ + ^(A(i,z)p(t,z))=0, (7) 
ot oz 

for the associated pdf, p{t, z). This formulation is similar to a pure convection problem (first-order 
partial derivatives in sample space). 
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As mentioned in the previous section, most of the time, the number of degrees of freedom (the dimen- 
sion of the state vector) is huge and one has to resort to a reduced (or contracted) description 1 13|. Con- 
sequently, one-particle pdf, p{t, z') (for a particle i) or two-particle pdf, p{t, z\ z^), etc, can be considered 
instead of choosing the N-particle pdf, p{t, z). Let be a reduced state vector (typically a one-particle 
state vector = {Zi, Z2, ■ ■ ■ , Zp) corresponding to the p variables attached to each particle). Then, the 
time evolution equations, in physical space, for this sub-system have the form: 

^=A(t,Z^,Y)), (8) 

where there is a dependence on the external variable Y (related to the particles not contained in Z^). In 
sample space, the marginal pdf p'^{t, z^) verifies 

^P-^^ + A-[{A\z^)p^{t,z^)]=0, (9) 
where the conditional expectation is defined by 



(A|z'-)= I A{t,z-,y)p{y\t,z-)dy=-^J A(t, z^ y) z^ y)dy. 



(10) 



Eq. (|9]l is now unclosed, showing that a reduced description of a system implies a loss of information and 
the necessity to introduce a model. 

In the present paper, and for reasons presented in the next section, the reduced system will be modelled 
by stochastic diffusion processes ll7ll9l [T2l[T4l . For such stochastic processes, the time-evolution equations 
for the trajectories of the process are SDEs written as: 

dZp = Ai{t, Z^{t)) dt + B,j{t, Z^{t))dWj, (11) 

where Wj = (Wi, . . . , Wn) is a set of independent Wiener processes [il4l and n is the dimension of 
the reduced state vector In Eq. (fTTT i. A — (Ai) is called the drift vector and B = (Bij) the diffusion 
matrix. SDEs require a strict mathematical definition of the stochastic integral II12II14I which is defined 
here in the Ito sense (these equations are referred to as Langevin equations in the physical literature). The 
corresponding equation in sample space for p^ {t, z'") is the Fokker-Planck equation 

¥- = -i^\AJt,z'-)p'-] + l-^^\D,,(t,z'-)p'-]. (12) 

where Dij = BuBij = {BB*)ij is a positive-definite matrix. In a weak sense (when one is only interested 
in statistics of the process), one can speak of an equivalence between SDEs and Fokker-Planck equations. 



2.2 Dimension of the state vector 

The dimension of the reduced state vector, Z^ , that is the number of particles N and the number of attached 
variables p for each particle, have to be determined (hereafter the upper-scripts r and R are dropped for the 
sake of simplicity). The first choice for N is done in line with current state-of-the-art models for single- 
phase flows Q. Indeed, when the particle relaxation time Tp is small, particles behave as fluid particles. 
In single-phase turbulence [15 |, only one-particle PDF models are sufficiently general to be applicable to 
complex flows. For this reason, our first choice is to retain a one-particle pdf description for the particle 
phase in the two-phase flows under consideration here {N ~ 1). 

The second choice is to select the specific variables attached to the solid particles. Again, a closer look 
at single-phase PDF models [TJ might be helpful. In single-phase flows at high Reynolds numbers, Kol- 
mogorov theory |fT6l tells us that, for a reference time scale dt in the inertial range, Lagrangian increments 
of the fluid velocity are well correlated whereas increments of the fluid acceleration are nearly uncorrected. 
This indicates that for dt belonging to the inertial range, the fluid velocity is a slow variable and the fluid 
acceleration is a fast variable which can be eliminated (fast variable elimination) ifTTl . Therefore, the state 
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vector should include position and velocity, i.e. (x/, U/) (p — 2). This is the starting point for Langevin 
equation models for fluid particle velocities [15"181. The model takes the form of a diffusion process with 
a drift term hnear in the velocity of the fluid seen IJSJ 



dxf^i — Uf,idt, (13) 
1 d{P) 



dUf^, = dt + G,, ({//,, - (Uf^,)) dt + VCo(e) dW,, (14) 

where {P) is the mean pressure field, (e) is the mean dissipation rate and Co is a constant given by Kol- 
mogorov theory (Co = 2.1). Gij is a matrix which depends on mean quantities, 

G,, = -^% + G^,. (15) 

where G°j is an anisotropy matrix (depending on mean quantities) and Tl stands for a timescale given by 
(k is the turbulent kinetic energy) 

Tl - (16) 

In the two-phase flow case, a similar reasoning fOl suggests to include the velocity of the fluid seen in 
the state vector that becomes (the fluid acceleration seen is a fast variable) 

Z = (Xp,Uj„U,). (17) 

This is different from the choice made in analogy to Boltzmann equation, when one considers only Z — 
(xp, Up) as in kinetic models I119II20I . Yet, we are dealing with particles being agitated by an underlying 
turbulent fluid and a (slow) variable related to the fluid, namely the velocity of the fluid seen, is explicitely 
kept in the state vector. With the kinetic choice, not only the derivatives of the fluid velocity seen have to 
be modelled but also the fluid velocity seen itself. 



3 Modelling turbulent dispersion 

With the present choice of the state vector, the stochastic process used to describe the system has been 
chosen, i.e. Z = (xp, Up, U^). Following the trajectory point of view mentioned in Section ITTl a time- 
evolution equation for U^ has to be proposed. This equation, together with Eqs ([T]i, will give the complete 
system of SDEs for the components of Z. Contrary to most Lagrangian models, which are often built in a 
discrete setting, the current model is written in continuous time, as Eq. (fTTT i. in order to be consistent with 
the proposed mathematical framework. 

From the physical point of view, a time-evolution equation for U^ amounts to modelling turbulent 
dispersion, an issue which is more complicated than turbulent diffusion. Indeed, particle inertia (Tp) and 
the effect of an external force field induce a separation of the fluid element and of the discrete particle 
initially located at the same point, as represented in Fig.[T] In the asymptotic limit of small particle inertia, 
Tp — ?• 0, and in absence of external forces, this separation effect disappears and the problem of modelling 
diffusion is retrieved, for which the stochastic model given by Eq. (fT4] l can be applied. For that reason, 
dispersion models (simulation of U^) are extensions of diffusion models (simulation of U/). 

An extensive description of the physical aspects of turbulent dispersion has been proposed elsewhere 
|I9]|2T|, so that only the key points used to derive the stochastic model are recalled in the next section. It 
is proposed to consider separately the physical effects of particle inertia and external forces. Two non- 
dimensional numbers have been introduced for that purpose: particle inertia is measured by the Stoke 
number St — t^/Tl, and external forces by ^ = \Ur\/u' , u' being a characteristic fluid turbulent velocity 
{u' — \/2kJi). The influence of these two effects on the characteristics of Us are : 

(i) in the absence of external forces = 0), only particle inertia plays a role. The characteristic, or 
integral, timescale of the velocity of the fluid seen, say T£ = 0) is expected to vary between the 
fluid Lagrangian timescale, Tl, in the limit of low St numbers, and the Eulerian timescale, Tg, in 
the limit of high St numbers. 
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(ii) Leaving out particle inertia, external forces creates mean drifts ^ 0) and induce a decorrelation 
of the velocity of the fluid seen with respect to the velocity of fluid particles. This effect is called the 
crossing trajectory effect (CTE) and is related to a mean relative velocity between particles and the 
fluid rather than an instantaneous one. 

In the model developed in the present paper, it is assumed that Te remains of the same order of mag- 
nitude as Tl, which seems actually a reasonable choice since there is little information for complex flows. 
Detailed models have been proposed for the effect of particle inertia [2X1^ but in the following it will be 
neglected, that is r£(^ = 0) = Tl. The representative picture is now sketched in Fig. |2] where only the 
mean drift induces separation. 

3.1 Langevin equation model 

Using the physical description of the CTE effect as due to a mean-drift (Fig. |2]i, Kolmogorov theory can 
be applied, as in the single-phase case, to suggest a dispersion model. Indeed, let us introduce v(t, r) = 
u/(io + T, Xo + w(ioj Xo)t + r) — u/(to, xq), the fluid velocity field relative to the velocity of the fluid 
particle F at time i„. Fig. |2] that is with u/(to, xq) = Us(to), then one can write that 

dU, = v(dt, (U^)di), (18) 

where (U^) = (Up) — (U^) is the mean relative velocity between the discrete particle and the surrounding 
fluid element. Then, the differential change, and so the Eulerian statistics, of the fluid velocity seen depend 
on the key variables of Kolmogorov (as the fluid velocity), that is < e > and v, and on the mean drift due 
to the CTE effect, but not on the instantaneous particle or fluid velocities. Since it is the mean velocity 
that appears in Eq. dTsl l, the Kolmogorov theory can then be applied 1 16|, to show that for high-Reynolds 
number flows and for a time increment dt that belongs to the inertial range, we have 

{dU,^,dUs^j)=D,j{dt), (19) 

where the matrix Dij is determined by the two scalars functions D\ | and D±^ through 

= D^5,, + [D\\~ D^_] nr„ (20) 

the separation vector r being in the direction of the mean relative velocity, r = (Ur ) / 1 (11^ ) | . The functions 
_D|| and are the longitudinal and transverse velocity correlation, respectively. Dimensional analysis 
yields that in the inertial range, one can write 

Dm . ie)dta, (^) , D.idt) = ie)dta, (^) . (21) 

For the two functions a|| and a±, there is no exact prediction, but in two limit cases they can be explicitely 
computed. On one hand, when the mean relative velocity is small, |(Ur)| ^ {{e)dt)^^^, for a given time 
interval dt, the statistics of the velocity of the fluid seen are expected to be close to the fluid ones, and thus 
a|| ~ a± ~ Cq. On the other hand, when the mean relative velocity is large, (|(Ur)| ^ {{e)dty^^), one 
can resort to the frozen turbulence hypothesis, and in that case (C is a constant) 

D^^idt) ~ C((e) (U,)di)2/3, D^idt) ~ ^C((e) {lJr)dt)^/\ (22) 

which shows that, in that limit, the two functions a\\{x) and aj_ (x) vary as a;^/"^. Then, the Langevin model 
is not supported as in the fluid case, since it will always give a velocity correlation linear in time for each 
components of D. Nevertheless, a useful approximation can be proposed. Indeed, if we freeze the values 
of the functions a|| and a± for a certain value of the time interval, say Atr, and write 

I (dt) - (^) -1 1 ) ^ idt) ^ (e) ) , (23) 
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a linear variation of I?| | {dt) and D± (dt), with respect to the time interval dt, is now obtained. The reference 
time lag may be the Lagrangian timescale which is the timescale over which fluid velocities are correlated. 
And since (e) Tl — k, we have 

D||(d<)^(e)dta|l(^^^^), D^(dt)^(e)dta^ (^Kl^^ . (24) 

This result suggests now a Langevin equation model which consists in simulating as a diffusion process. 
As explained above, this model is only an approximate model having less support than in the fluid case. 
Indeed, the Langevin model does not yield the correct spectrum (in the limit of large relative velocity 
or frozen turbulence). However, for engineering purposes, where the macroscopic behaviour is the real 
subject of interest, the important properties are the integral time scales rather than the precise form of the 
spectrum. Thus, Langevin models are "reasonable compromises" between simplicity and physical accuracy 
at the moment. It is also clear that much work remains to be done to improve stochastic models. 

It can be shown [9 1 that the general stochastic differential equations for the fluid velocity seen process 
have the form (X stands for fluid fields) 

dUs,, ^ A,{t,Z,{Z),{X))dt + B,,it,Z,{Z),{X))dWj, (25) 

where the drift vector, A^, and the diffusion matrix, B^, have the form 

-7^{U,j-{U,,,))dt 



+ y (e) {Cohk/k + ^{hk/k - 1)^ dW,. (26) 

The CTE has been modelled by changing the timescales in drift and diffusion terms according to Csanady's 
analysis. Assuming, for the sake of simplicity, that the mean drift is aligned with the first coordinate axis 
(the general case is discussed elsewhere ||9]), the modelled expressions for the timescales are, in the 
longitudinal direction 

T*{£ — 0) 

T£ , = , ' =, (27) 




and in the transversal directions (axes labelled 2 and 3) 



T*(f = 0) 

T£,2 - r£ 3 = , ' = . (28) 




2fc/3 



In these equations f3 is the ratio of the Lagrangian and the Eulerian timescales of the fluid, /S — Tl /Te, and 
r£(^ = 0) represents the Lagrangian time-scale in the absence of mean drifts but accounting for particle 
inertia. As mentioned at the end of the previous section, particle inertia effect are neglected in the present 
work f|211 and we therefore assume that r£(^ = 0) = T^. In the diffusion matrix, a new kinetic energy 
has been introduced (bi — TL/Tj^^i) 

In the absence of mean drifts, the stochastic model for reverts to the Langevin equation model used 
in single-phase PDF modelling |7| and is thus free of any spurious drift by construction. Finally, it must 
be emphasised that the derivation of a satisfactory model (that is respecting a number of well-established 
constraints) for particle dispersion remains an open issue. 
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3.2 Modelling two-way coupling 

In order to account for the influence of the particles on the fluid, a new term is added in the momentum 
equation of the fluid velocity, see Eqs (|4|l, and the fluid velocity seen, 

dUs,, = [As4t, Z, (Z), (X)) + Ap^sAt, Z, (Z))] dt + Bs,^,{t, Z, (Z), (X)) dW,. (30) 

The exact expression for this acceleration, Ap^s,i(i, Z, (Z)), which is induced by the presence of the 
discrete particles, is not a priori known. The underlying force corresponds to the exchange of momentum 
between the fluid and the particles, but should not be confused with the total force acting on particles since 
the latter includes external forces such as gravity. The effect of particles on fluid properties is expressed 
directly in the stochastic equation of with a simple stochastic model. The force exerted by one particle 
on the fluid corresponds to the drag force written here as 

^p^j = -mp i^, (31) 

Tp 

where nip is the mass of a particle. The total force acting on the fluid element surrounding a discrete particle 
is then obtained as the sum of all elementary forces, Yp^j, and the resulting acceleration is modelled here 
as ID 

. _ CKpPp Up,i — Us,i 

afPf Tp 

Eq. (l30l l is justified by the assumption that the mean transfer rate of energy and energy dissipation (e) 
is changed by the presence of particles, but the nature and structure of turbulence remains the same. There- 
fore, Eq. (l30l l is written by adding an acceleration term, Ap^s, to account for the presence of particles, 
while the same closures as in the one-way coupling case will be used for the drift vectors and the diffusion 
matrices, where, once again, the mean fields (e), {U'j), . . . are modified by the presence of the particles. 
Indeed, the drift vectors and the diffusion matrices not being affected by the nature of turbulence, remain 
unchanged. In opposition to the previous hypotheses, recent results of direct numerical simulations in the 
field of turbulence modulation by particles (in isotropic turbulence) |4| seem to indicate that there is a 
non-uniform distortion of the energy spectrum. This could mean that, contrary to our previous assumption, 
the nature and structure of the energy transfer mechanisms of turbulence are modified by the presence of 
particles. There is no precise 'geometrical' knowledge on the structure of turbulence in the presence of 
discrete particles and this makes it extremely difficult to isolate the important variables in order to modify 
the theory of Kolmogorov (which is used in our closures). This problem is out of the scope of the present 
paper and it remains an open question. Then, the final set of equations for the velocity of the fluid seen are: 

dUs^. = ~^^dt+iiUp.,)-{Uj,))^dt~^^^I^^^^dt 
Pf dxi dxj afpf Tp 

^ - {Us.,- {Uf,,)) dt 



+ J(e) {c^hk/k + '^{hk/k - l)^ dW,. (33) 

It is seen that the resulting Langevin equation, which is believed to represent the simplest model for two- 
phase flows, contains a diagonal but non-isotropic diffusion matrix, Bs.ij — Bs^i&ij- It is also worth 
emphasising that the closure relations put forward just above reflect modelling choices. For instance, in 
the two-phase flow case, the isotropic form of the diffusion matrix cannot be obtained anymore, but it is 
chosen to select among different possibilities, a diagonal diffusion matrix. 
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3.3 Equivalence with the PDF approach 



According to the arguments developed in Section |2n the complete set of SDEs (for the state vector Z = 
(xp,U/,U,)), 

dxp^i = Up^i dt (34a) 
dUp.^ = — {Us,^ - Up^i) dt + g, dt (34b) 

Tp 

dUs,^ = Z, (Z), (X)) + Ap^sAt, Z, (Z))] dt + Z, (Z), (X)) dW,. (34c) 

is equivalent to a Fokker-Planck equation given in closed form for the corresponding pdf p{t\ y-p, Vp, V^) 
which is, in sample space 

9p ... dp _ __d 

^ ([A... + {Ap^sAyp^^P^s)] P) + Lot/^t. {[BsBJ],,p). (35) 



The equation for the Eulerian pdf and the resulting mean-field equations can be found in ||8l . 

4 Numerical Issues 

The theoretical model developed in Section [3]represents a PDF model for the particle phase only. It does 
not contain any description of the continuous phase. It is possible to extend the PDF description to both 
the fluid and particle phases [[8J, which may be useful for theoretical and consistency analysis. However, 
at the moment, this complete PDF approach is limited for practical calculations, and, in the present work, 
a classical second-moment approach is followed for the continuous phase. The complete numerical model 
is therefore an hybrid method and corresponds to a classical approach referred to as Eulerian/Lagrangian 
in the literature, as mentioned in Section[T] As one can see from Section|3] the terminology is not actually 
adequate to describe the complete model (it would be better to talk of a Moment/PDF hybrid approach), but 
corresponds to the numerical approach. Indeed, from the numerical point of view, the fluid phase is mod- 
elled by mean fields, obtained by solving partial differential equations on a grid with an Eulerian approach, 
while the particle phase is modelled by a large number of Lagrangian particles distributed in the domain 
and whose properties are obtained by solving stochastic differential equations. It is worth emphasising that 
these particles are now stochastic particles, or more precisely samples of the underlying pdf, rather than 
precise models of the actual particles. The overall numerical method is therefore an example of Monte 
Carlo particle-mesh techniques. 

The numerical (particle-mesh) approach involves many issues. Some of them have already been treated 
in classical textbooks |22|, but only for deterministic equations. The stochastic nature of the present 
equations brings in specific aspects and raises new questions. In that respect, the purpose of this section 
is not to give a comprehensive description of all issues. It is more to give an overview of the numerical 
method, pointing out important issues and, in particular, those that, in our opinion, require additional 
work. More precisely, issues that have not always been investigated or may have been overlooked (such as 
consistent discrete averages, Section l431 and mass-continuity constraint, Section l4~4] i. are developed more 
in detail. 



4.1 General Algorithm 

The flow-chart of the code is shown in Fig. |3] At each time step, the fluid mean fields are first computed 
by solving the corresponding partial differential equations (RSM model) with a classical finite volume 
approach. The Eulerian solver then provides the Lagrangian solver with the fluid mean fields that are 
necessary to advance particles properties. In the Lagrangian solver, the dispersed phase is represented by 
a large number of particles and, as proposed by the model, the state vector attached to each particle is 
Z — (xp,U/,Us). Once particle properties have been updated, and in the case of two-way coupling, 
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where particles modify the fluid flow, source terms accounting for momentum and energy exchange be- 
tween the two phases are also calculated and are fed back into the Eulerian solver for the next time step 
computation. It is then seen that the two solvers are only loosely-coupled. This may lead to numerical 
difficulties when the particle loading is increased, consequently the source terms become important and 
the system of equations stiff. However, our present aim is to model moderate particle loading phenomena, 
indeed particle-particle collisions have been neglected. In that range, particles can still modify the fluid 
flow in a noticeable way but source terms remain small enough so that the loosely-coupled algorithm can 
still be retained. 

As previously explained, the particle properties are modelled by a vectorial SDE written as 



where / is a general function depending on the model and X stands for fluid fields. It is worth emphasising 
that the drift and diffusion coefficients depend on statistics derived from the pdf that is implicitly calculated. 
Therefore, these SDKs are different from standard ones lfT4l|26l . Updating particles properties implies 
three steps: (i) projection of (/(Z)) and (X) at particle positions, (ii) time integration of Eq. (l36b . and 
(iii) averaging to compute the new values of (/(Z)) (for stationary flows, such as the one considered later 
on, ensemble averages computed in every cell are then averaged in time, once the stationary regime has 
been reached. This time-averaging procedure is very helpful to reduce statistical noise to a negligible level 
ll23H25j ). Since averaging is basically the reverse operation of projection L22J . these three steps correspond 
to two main issues: 

(i) the first concerns the derivation of accurate numerical schemes for the time integration of Eq. ( [36l l. 
Due to the non-linear nature of the equations, this is still a difficult point It26ll27l and, moreover, 
physical constraints should be respected [28 J . This issue is briefly developed below. 

(ii) The second issue is related to the exchange of information between the grid-based Eulerian variables, 
located at cell centres and particles which are continously distributed in the domain. At the moment, 
a NGP (nearest grid point) technique 1221 is used, this represents the simplest choice but also the 
best one in terms of spatial error ||23l. This is an important and attractive issue to investigate for 
particle-mesh methods with in the case of unstructured meshes and taking into account boundary 
conditions. 

4.2 Time-integration of SDEs 

Since we are interested in the numerical approximation of statistics derived from particles, a weak numer- 
ical scheme lIZTll (converging in law) is under consideration. A numerical scheme is said to be of order of 
convergence r in time, in the weak sense, if, for any sufficiently smooth function 



where C is a constant and Z represents the numerical approximation of Z. The numerical scheme 
used in the present calculations is detailed in ll28l . and consequently, in the present paper, only some 
points of particular importance are emphasised. Eq. ( l36l l must be understood in the ltd sense and it is 
fundamental that numerical schemes respect the Ito definition of the stochastic integral, in order to avoid 
any inconsistency problems ||29l . The weak numerical scheme is of order 2 in time, unconditionally stable 
but still explicit 1 28 1 . Another important issue is the numerical fulfilment of physical limits 1 9 1 . Indeed, in 
practical engineering calculations of complex flows, it may occur that, locally, one has Ai Tp, or even 
At ^ T*, Tp, that is the time-step becomes much larger than the characteristic time scales of the system, 
Eqs ( l34l) . In the first case, one should have that Up — > Ug and, in the second case, the model expresses a 
pure diffusive behaviour in space [9 , 28 1 



dZ, = Mt, Z, (/(Z)), {X))dt + B,,{t, Z, (/(Z)), {'K.))dWj 



(36) 



|(/(Z))-(/(z^*))|<c(A^)^ 



(37) 



dxp,i = {xf^i) dt + {Bs,ij Tl)dWi. 



(38) 



It is important that the numerical scheme is consistent with these continuous limits. 
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4.3 Discrete representation and numerical averages 

Since averages are fundamental in the construction of PDF models, it is useful to clarify the correspondence 
between the averages (defined as the mathematical expectations) and Monte Carlo estimations, which are 
used in the code. In polydispersed cases, even when pp is constant, the mass of each particle can be different 
because of different diameters. This suggests that even for constant density particles, the natural definition 
or understanding of a mean quantity is the mass-weighted average. This kind of choice is somewhat 
analogous to the Favre mean definition for compressible single-phase fluid flows. 

To justify this, we start by introducing a Lagrangian mass density function F^{t; y-p, Vp, 'i'p) where 

i^^(t;yp,Vp,*p)dypdVpd*p, (39) 

is the probable mass of discrete particles in an infinitesimal volume in sample space. As a matter of 
fact, attention is focused on Eulerian averages (at a point {t, Xp) fixed in time and space), for which the 
analogous Eulerian mass density function is defined by 

F^(t,x;Vp,t/.p) = F^(t;yp-x,Vp,V'p) 

f r (40) 

= / ^ (i;yp,Vp,Vp)(5(x-yp)dyp, 
where is normaUsed by {{pp) is the expected density) 

ap(i,x)(pp)(t,x) = J F''{t,^;Vp,tp,^)dYpd'iPp. (41) 

ap represents the probability to find particles at a given time and position, in any state. The Eulerian mass 
density function being defined, we can introduce a general average for a quantity H{XJ{t), 4>{t)) 

ap(<,x)(pp)(t,x)(ffp)(i,x) = y"i7p(Vp,*p)F^(i,x;Vp,*p)dVpdtAp. (42) 

The Lagrangian mass density function can be written from a discrete point of view as 

N 

Fk{t-,ypyp,'^p) = ^m*<5(yp - ^p{t)) ® 5(Vp - U;(t)) ® 5{xPp - (43) 

1=1 

where is the mass of the particle labelled i and N is the number of samples. From Eq. ( |39] l. the discrete 
Eulerian mass-density functions is 

1 ^ 

F§{t, Xp; Vp, V^p) = _ ^ m^<5(Vp - Uj,(t)) ® Si% - 0;(t)) , (44) 

^ i— 1 

6Vx being a small volume around point x. Then, a numerical approximation of Eq. (|4TI ) is 

apM {pp)^ ^J^^\ (45) 
and the numerical approximation of a particle mean quantity is 

{Hp) ~ Hp^N = — jv : • (4o) 

Convergence of the discrete approximation is ensured by the Central Limit Theorem which shows that 
there exists a constant C such that, when N +oo, 

{Hp^N)^{Hp) and {{{Hp) ^ Hp^ri?) < (47) 

It is therefore seen that the convergence of the underlying pdf is not in a strong sense but in a weak sense, 
or to be more precise in law I.14J, since it is in fact the mean value of functions of the stochastic process Z 
that converges as iV ^ +oo, 

Hn.p = {H{Zn)) > {H{Z)). (48) 
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4.4 Pressure correction 



It has been shown in Section [331 that there is a correspondence SDE - Fokker- Planck equation. From the 
pdf equation, mean particle fields can then be extracted HI. In other words, every particle stochastic model 
is consistent with a certain Eulerian model ||6l|9] as in single-phase PDF models |30|. 

With the definition of the mean particle velocity field given in the previous section, the corresponding 
particle continuity equation is (density is constant) 

d d 

^KPp) + -g^^iap Pp{Up^^)) = 0. (49) 

For each time step in the Lagrangian solver, the mean fields Up and {Up.i) are computed from particle 
location and velocity, Xp and Up, using the numerical approximations given in Eqs ( |45]) and (|46] |. Here, 
we propose to modify particle velocities (not locations) so as to enforce the mean continuity constraint, by 
adding a pressure-correction field as a potential 0. The corrected particle velocity field is then 

{Up.,) = {Up,^) - (50) 

where cf) is calculated from the Poisson equation 

9 , d4) . d , . d , ITT x-, 

d^^'^^P^d^^ ^ dl^'^^P'^ + a^KPp(^P.)) • (51) 

The mean velocity correction term is then applied to each particle velocity. 

This pressure-correction term used here for the particle velocity is of course similar to the classical 
pressure-correction step applied in the Eulerian solver for the fluid. Yet, it is often overlooked in La- 
grangian calculations. If we consider the complete algorithm, it is then seen that there are now two pressure- 
correction steps due to the two mean-continuity equations, one for the fluid and one for the particles. This 
is also a consequence of the loosely-coupled algorithm. 



5 Numerical investigation 
5.1 Experimental setup 

The experimental setup is typical for pulverised coal combustion where primary air and coal are injected 
in the centre and secondary air is introduced on the periphery. Fig. ID 

This is a typical bluff -body flow where the gas (air at ambient temperature, T = 293 K) is injected in 
the inner region and also in the outer region where the inlet velocity is high enough to create a recirculation 
zone downstream of the injection (two honeycombs were used in the experiment in order to stabilise the 
flow so that no swirl was present). Solid particles (glass particles of density pp = 2450 kg/m?) are then 
injected from the inner cylinder with a given mass flow rate and from there interact with the gas turbulence. 
This is a coupled turbulent two-phase flow since the particle mass loading at the inlet is high enough (22%) 
for the particles to modify the fluid mean velocities and kinetic energy. This is also a polydispersed flow 
where particle diameters vary according to a known distribution at the inlet, typically between dp — 20/im 
and dp ~ 110/xm around an average of dp ^ 60 jim. 

Experimental data are available for radial profiles (the flow is stationary and axi-symmetric) of different 
statistical quantities at five axial distances downstream of the injection {x = 0.08, 0.16, 0.24, 0.32 and 0.40 
m). These quantities include the mean axial and radial velocities as well as the fluctuating radial and 
axial velocities for both the fluid and the particle phase. Axial profiles along the axis of symmetry for 
these quantities have also been measured. All the data was gathered using PDA measurement techniques. 
Further details on the experimental setup and the measurement techniques can be found in Ishima et al. 

USD. 

The 'Hercule' experimental setup is a very interesting test case for two-phase flow modelling and 
numerical simulations where most of the different aspects of two-phase flows are present. The particles are 
dispersed by the turbulent flow but in return modify this one. Furthermore, the existence of a recirculation 
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zone where particles interact with negative axial fluid velocities constitutes a much more stringent test case 
compared to cases where the fluid and the particle mean velocities are of the same sign (the problem is then 
mostly confined to radial dispersion issues). 



5.2 Results and discussion 

All the results were obtained by using the ESTET 3.4 software on a HP-C3000 workstation. In all numerical 
computations, the axi-symmetry property was used: a two-dimensionnal curvilinear mesh with 74 x 3 x 
142 nodes was generated. The sensitivity to the various parameters of the numerical investigation was 
accurately studied. In particular, independence with respect to the time step was checked. A uniform time 
step. At = 10^^ s, was then used in all computations. 

The computations were carried out with a Rij — e turbulence model, which is based on the standard 
1PM model II32II33I . Actually, this choice is satisfying from the point of view of the consistency with the 
stochastic model. It is known that there is a rigorous correspondence between the Lagrangian stochastic 
models and the second-order closures in the case of turbulent single-phase flows 1301 . 

In the two-phase flow calculations, particles were injected when the single-phase flow stationary regime 
was reached (as the limit of the unstationary regime) before the introduction of the discrete particles in 
the domain. About 1000 time steps were computed for the single-phase problem. Around 400 to 500 
additional time steps were needed to reach the stationary regime for the two-phase flow situation (around 
14000 particles were at this stage present in the domain). Statistics extracted from the particle data set were 
then averaged in time (for about 1000 time steps) to reduce the statistical noise. 

The computational performances are shown in Table[T] Normally, Lagrangian algorithms require much 
more computational time than the Eulerian eddy-viscosity models |9|. In this case, for the same num- 
ber of computational elements (either mesh points or nodes), they appear comparable. The computational 
requirements for the Eulerian solver is increased due to the use of a full second-order turbulence model 
which implies the numerical resolution of 6 coupled partial differential equations for the fluctuating veloc- 
ities (added to the 3 equations for the mean momentum) compared to only 1 for eddy-viscosity models. 

The experimental set of measures provides data both along the axis and in cross sections at various 
points in the domain. The comparison is made in all directions and at all cross sections of measures. The 
cross section at a; = 0.16 is located within the recirculation zone while the cross section at a; = 0.4 is 
located downstream of the limit of the recirculation zone. 

The overall agreement between experimental data and the computed profiles is good. In particular, the 
particle fluctuating velocity is well reproduced both in shape and in magnitude. 

In Fig. |5]and|6] the mean fluid and particle velocities along the axis are shown. It is noticeable that the 
comparison between the computed results and the experimental findings for the two-phase flow in presence 
of two-way coupling is worse than in single-phase computation. The same effect characterises both the 
mean fluid and the particle profiles. They results less well reproduced in two domain zones, although the 
qualitatively agreement remains good. The point of recirculation is overestimated and the velocity slope 
after it is underestimated. This effects indicates the necessity of further studies on the coupling between 
particles and the fluid. It is worth noting that these effects are limited to the behaviour along the axis. In 



Fig. [SifTOl first two statistical moments of the particles velocity ((Up), (Wp), y (u'p)) are shown, without 

smoothing. The difference between experimental data and computed results at the axis (x — 0) does not 
influence the computation in the rest of the domain. In Fig. |7]we show the profiles of the fluid mean 
axial velocity, where an analogous behaviour is present, with a satisfactory agreement in the whole domain 
except the values on the axis. 

6 Conclusions 

In this paper, a theoretical and numerical model for particle turbulent polydispersed two-phase flows has 
been presented. The theoretical model is a PDF model and, in practice, appears as a Lagrangian stochastic 
model. It consists in the simulation of a large number of stochastic particles which simulate the behaviour 
of real particles dispersed in the fluid. Each particle is defined by a set of variables and the selection 
of these state variables represents an important choice from the physical point of view. At present, the 
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state variables attached to each particle include particle position, particle velocity and the fluid velocity 
seen. The present model is developed as a diffusion stochastic process for the velocity of the fluid seen. 
This is similar to single-phase turbulence, but the extension to the two-phase flow case requires additional 
assumptions in the application of the Kolomogorov hypotheses, as detailed in Section [3] A specific point 
is that, in the case of two-way coupling, an extra term is needed in the stochastic equation for the velocity 
of the fluid seen in order to be consistent with the mean field equations for the fluid phase. 

From the numerical point of view, an hybrid Eulerian/Lagrangian, or moment/Monte Carlo, approach 
is discussed. At each time step, the fluid phase is computed with an Eulerian code which provides the 
Lagrangian module with mean fluid quantities. The particles are then tracked and source terms representing 
the momentum and kinetic energy exchanges are evaluated to be included in the Reynolds stress equations. 
This corresponds to a classical approach, but new aspects have been emphasized. In particular, apart from 
considerations on numerical schemes and the evaluation of particle means, the necessity of a correction to 
satisfy particle continuity equation has been stressed. 

The interests and capabilities of the model have been illustrated by the computation of a test case repre- 
sentative of an engineering situation. Numerical predictions are in good agreement with the experimental 
ones and can be regarded as a validation of the model. 

Some of the current developments to this work aim at improving numerical aspects (variance reduction 
technique for the computational efficiency, new methods to compute statistical averages) and at improving 
the physics of the model in the near-wall region (boundary layer). 
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CPU time for 1000 nodes/time step 


CPU time for 1000 particles/time step 


Eulerian solver 


0.20 s 




Lagrangian solver 




0.17 s 



Table 1 : Computational performances 
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fluid particle trajectory 
Figure 1 : Fluid element and particle paths 
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Figure 2: Mean fluid and parlicle paths 
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Figure 3: Sketch of the algorithm for one time step 



17 



Uj = 4 m/s 
Re=4500 



x=0.08 m 



x=0.16m 



x=0.24 m 



x=0.40 m 




L=l m 



Figure 4: The 'Hercule' experimental setup. The mean streamlines are shown for the fluid (solid lines) 
and the particles (dashed lines). Two stagnation points in the fluid flow can be observed (Si and 82)- 
Experimental data is available for radial profiles of different statistical quantities at five axial distances 
downstream of the injection (x = 0.08, 0.16, 0.24, 0.32, 0.40 m) (experimental data is also available on the 
symmetry axis). 
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Figure 8: Profiles of mean axial particle velocity 
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